Quantum Monte Carlo Analysis of Exchange and Correlation in the Strongly 

Inhomogeneous Electron Gas 
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We use variational quantum Monte Carlo to calculate the density-functional exchange-correlation hole Uxc, the 
exchange-correlation energy density e^c, and the total exchange-correlation energy Exc, of several electron gas 
systems in which strong density inhomogeneities are induced by a cosine-wave potential. We compare our results 
with the local density approximation and the generalized gradient approximation. It is found that the nonlocal 
contributions to exc contain an energetically significant component, the magnitude, shape, and sign of which are 
controlled by the Laplacian of the electron density. 
PACS: 71. 15. Mb, 71.10.-w, 71.45. Gm 
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Kohn-Sham density-functional theory (DFT) shows 
that it is possible to calculate the ground-state proper- 
ties of interacting many-electron systems by solving only 
one-electron Schrodinger-like equations. The results are 
exact in principle, but in practice it is necessary to ap- 
proximate the unknown exchange and correlation (XC) 
energy functional, i?a;c[n(r)], which expresses the many- 
body effects in terms of the electron density n(r). The 
current popularity of density-functional methods in con- 
densed matter physics, quantum chemistry, and materi- 
als science reflects the remarkable success of fairly simple 
approximate XC energy functionals. 

In the local density approximation (LDA) , the XC hole 
n^c about an electron at r is assumed to be the same 
as in a uniform electron gas of fixed density n = n(r). 
This works surprisingly well, even though real solids and 
molecules have strongly inhomogeneous electron densi- 
ties. Almost all attempts to go beyond the LDA have 
been based on studies of XC in the weakly inhomoge- 
neous electron gas. The widely used generalized gradi- 
ent approximation (GGA) |^-^ incorporates information 
about first-order gradient corrections to the LDA, as ob- 
tained in the limit of slowly varying electron densities. 
This is done in a controlled way such that, for exam- 
ple, Tixc satisfies known sum rules but without any 
guidance from the behavior of n^c in the strongly inho- 
mogeneous regime. This may explain why current GGAs, 
although better than the LDA, are not consistently able 
to deliver the very high accuracy of ~0.1 eV required 
to study many chemical reactions. To provide a new 
direction in the search for more accurate and reliable ap- 
proximate functionals, we have used variational quantum 
Monte Carlo (VMC) ||] to explore exchange and corre- 
lation in the strongly inhomogeneous electron gas. 

In this letter we calculate the central quantities in 
DFT: the XC hole rixc, the XC energy density Cxc, and 
the total XC energy E^c- We study an electron gas of 
average density nP corresponding to rs=2 and induce 
strong inhomogeneities by applying a cosine- wave poten- 



tial yqcos(q.r). Fixing Vq at 2.08ep, we consider several 
values of g < 2.17fc^, where and kp are the Fermi 
energy and wave vector of a uniform electron gas with 
r,=2. 

The electron density in our systems varies strongly and 
rapidly on the scale of the inverse local Fermi wavevector 
/ci?(r)~^ = (37r^n(r))~-'^/^, producing a strikingly nonlo- 
cal behavior of Uxc that cannot be described by semilo- 
cal corrections to the LDA. We show, however, that the 
resulting LDA errors in Cxc have a dominant and ener- 
getically significant component, the magnitude, shape, 
and sign of which are controlled by V^n(r). The GGA is 
unable to correct the LDA errors in E^c resulting from 
this component in an adequate way, a deficiency that 
severely restricts its accuracy. The relevance of Lapla- 
cian terms has been pointed out previously [^[js], but 
our calculations provide the first quantitative evidence 
of their importance in strongly inhomogeneous systems 
and for Cxc- 

The starting point of our calculations is the adiabatic 
connection formula In Hartree atomic units (e = 
m = 47reo = h = 1), this expresses i?a;c[rt] as the volume 
integral of an XC energy density Cxc defined by: 



dr 



, n{r)nxc{r,r') 



(1) 



where Uxc is the coupling-constant-averaged XC hole. 
The definition of Cxc used in the construction of the GGA 
differs from ours by an integration by parts that does not 
affect the integrated Exc- We favor Eq. (|l|), however, be- 
cause of its clear physical interpretation and because it 
aids comparison with the LDA XC energy density, which 
is constructed from an approximate hole. 

The XC hole Uxc is obtained via a coupling-constant 
integration, 

n(r)n(r') -I- n(r)n^c(r, r') 

= /' d\ (*a| E E -^(^ - - , (2) 
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where ^\ is the antisymmetric ground state of the Hamil- 
tonian H'^ = T + XVee + associated with couphng 
constant A. Here T and Vee are the operators for the 
kinetic and electron-electron interaction energies, and 
= T,V^{ri) is the one-electron potential needed to 
hold the electron density n^{r) associated with equal 
to n(r) for all values of A between and 1. At full cou- 
pling (A=l), y^(r) coincides with the external potential 
of the system. For other A values V''^ must be determined. 

Our VMC method for calculating n^c and Cxc from 
Eqs. dH-p ) is a generalization of the scheme used by Hood 
et al. It amounts to treating both '^^ and vari- 
ationally and determining the variational parameters by 
simultaneously minimizing the variance of the local en- 
ergy [[l2yi3| ] and the deviation of n'*' from n This 
is achieved by adding a term proportional to |n^-n| to 
the variance of the energy and minimizing the resulting 
penalty function, which attains its absolute minimum of 
zero only when satisfies the Schrodinger equation as- 



sociated with A and n{r) 



The VMC calculations are performed for a finite spin- 
unpolarized electron gas in a FCC simulation cell sub- 
ject to periodic boundary conditions. We generate the 
ground-state density by applying the cosine- wave poten- 
tial to the A=0 system and solving the self-consistent 
Kohn-Sham equations within the LDA to obtain the 
single-particle orbitals We then define the result- 
ing density to be the ground-state density of the fully 
interacting system We consider potentials with 

q = l.llkp, 1.55kp, and 2.17A;^, and cells containing 64, 
78, and 69 electrons, respectively. The large amplitude of 
the cosine-wave potential ensures that these systems are 
in the strongly inhomogeneous regime, |n(r)— n"|/n" « 1. 

The adiabatic calculations were performed using 6 
equidistant values of A in the range [0, 1], and with the 
following Slater- Jastrow ansatz for "i^: 



cxp 



i>j 



, (3) 



where = jr^ — rj\ and and are spin-up and 
spin-down Slater determinants constructed using the ex- 
act Kohn-Sham orbitals 4>i. The two-body term u^. ^. 
contained both a fixed part and a variational part as dis- 
cussed in The one-body term the potential V^, 
and the electron density were all expanded in plane 
waves. At each A we used a total of 20 variational pa- 
rameters in and and up to 7 coefhcients in the 
plane- wave expansions of and v^. The optimization 
of the parameters in '^■^ and was performed using 
96000 statistically uncorrelated electron configurations. 
This was sufficient to reduce the root mean square de- 
viation of n^(r) from n(r) to less than 0.5% of n(r) for 
all values of A and all systems. After optimization, ex- 
pectation values were calculated using the Monte Carlo 



Metropolis method with 10^ independent configura- 
tions of all electrons. Throughout we used the modified 
electron-electron interaction described in , which vir- 
tually eliminates the finite-size errors arising from the 
long range of the Coulomb potential. The statistical 
errors were negligible except in the tails of n^c in low- 
density regions, where they were less than 3%; this is 
much smaller than the differences between the VMC and 
LDA XC holes. 

The largest systematic errors are caused by the finite 
size of the system and the approximate nature of 4*^. 
These errors combine such that, even in a homogeneous 
electron gas, e^j^*-^ ^ ^x!?^- To circumvent this prob- 
lem we performed additional VMC calculations for finite 
homogeneous electron gases with N—QA and = 0.8, 
1, 2, 3, 4, 5, 8, and 10. The results enabled us to con- 
struct a Perdew-Zunger (PZ) parameterization of the 
VMC XC energy per electron of a finite uniform electron 
gas with A^=64. This parameterization was used to cal- 
culate e^c^ and e^c^ in all systems studied, ensuring 
that these quantities were exactly equal to e^/^*^ in any 
finite homogeneous system with Af=64. This procedure 
largely eliminates the systematic errors in the calculated 
differences between e^^'^'^ and e^f'^ reported below. 

In Fig. |l] we show snapshots of the deformation of 
"afc*^'^ around an electron moving in the q=\.bbkF sys- 
tem along a line parallel to q, the direction of maximum 
inhomogeneity, from a density maximum towards the tail 
of nir). The XC hole is plotted as a function of r' around 
a fixed electron at r, with r' ranging in a plane parallel to 
q. Also shown is the corresponding LDA hole n^/r'^ 
At the density maximum (not shown) both n^^'^^ and 
^xF^ are centered around the electron. However, un- 
like n^^^, which is always spherically symmetric, n^c^'^ 
is contracted in the direction of inhomogeneity and less 
compact. As the electron moves away from the density 
maximum to a point on the slope (top panel), the non- 
local nature of n^^'^^ becomes manifest. While n^^^ 
is still centered around the electron and is rather dif- 



fuse, 



^VMC 



lags behind near the density maximum and is 



much more compact. The nonlocal behavior of n^^^'~^ be- 
comes remarkable at the density minimum. Here n^^^''^ 
has two large nonlocal minima, each centered around a 
density maximum ~2.80 a.u. away from the electron. 
The LDA hole at this point is much more long ranged 
than n^.^'^'-^ and is spread over the whole system in order 
to satisfy the sum rule / (ir'n^^'^(r, r') = —1. 

This striking nonlocality of n^^'^'^ also occurs in the 
other two systems we considered. Clearly, semilocal cor- 
rections are unable to significantly improve the LDA de- 
scription of the XC hole in our systems, and fully non- 
local approximations are required. We found, however, 
that despite the strong nonlocality of n^^^'-^ , the LDA 



^VMC 

-'XC 



quantity, the Laplacian of the electron density. 

In Fig. H we show e^^^ — e^^^ for two of the strongly 
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inhomogeneous systems studied, where e^f'^ is calcu- 
lated using the exact ground-state density n(r). The re- 
sults are plotted along a line parallel to q (we call this 
direction y). Also shown are n(r) and V^n(r), plotted 
along the same line. It is apparent that the shape, magni- 
tude, and sign of the LDA errors in Cxc closely follow the 
shape, magnitude, and sign of V^n(r). The LDA errors 
in Cxc are large and negative in regions where V^n(r) is 
large and negative (around density maxima), and large 
and positive in regions where V^n(r) is large and pos- 
itive. The GGA XC energy density is not defined via 
Eq. (|l|) and so GGA errors are not shown. 

The VMC values of the integrated E^c are shown in 
Table |, along with the differences AE^^^ = E^^^ - 
E^^^^ and AE^^"^ = Eg'^^-E^J^^ (The version of the 
GGA used here is due to Perdew, Burke, and Ernzerhof 
[Ql ) . The LDA errors in E^c reflect the profound effect of 
the Laplacian errors in Cxc and change sign from positive 
(for the system ) to negative (for the two other 

systems) as q increases and the negative contributions to 
Acxc, which occur where V^n(r) < 0, become dominant. 
The GGA corrections are by construction always negative 
they improve E^^^ for the q—l.llkp system but 
worsen it for the two other systems. 

Since any smooth density n(r) may be expanded as 
a power series, the XC energy density functional may 
always be written as 

exc{r, [n]) = exc{r, n(r), Vin(r), ViVjn(r), . . .) . (4) 

The GGA may be viewed as an attempt to approximate 
the right-hand side of Eq. (|4|) as a nonlinear function of 
n(r) and Viri(r). In the case of the exchange energy, a 
uniform scaling argument (l7[ shows that 



Fx{p,l,...)et 



LDA 



(5) 



where Fx is an enhancement factor, Cx^^ is the LDA ex- 
change energy density, and p — \Vn\/ {2kF{'r)n{r)) and 
I — V^n/(4fc|n(r)n(r)) are a dimensionless gradient and 
a dimensionless Laplacian, respectively. We have cal- 
culated Fx exactly for our systems and investigated its 
dependence on p and /. Figs. ^ and are scatter plots 
showing values of Fx plotted against values of p and 
respectively. 

The two-valued nature of Fig. ^ arises because the 
mapping from density gradient to position is two valued 
in our systems. For example, the density gradient is zero 
both at the minimum and the maximum of the density, 
where the required corrections to the LDA XC hole are 
completely different, as can be inferred from Fig. |l|. In 
sharp contrast, Fx is a simple and almost unique func- 
tion of I. This suggests that the inclusion of Laplacian 
terms may allow the construction of simpler and more 
accurate approximate functionals. In most gradient ex- 
pansions the Laplacian terms allowed by symmetry are 
transformed into |Vn(r)|^ terms via an integration by 



parts. This is only possible when the dependence on I 
is linear, however, which is not the case here as can be 
seen from Fig. |3|3. Furthermore, the integration by parts 
destroys the physical interpretation in terms of the XC 
hole and so hinders further progress. 

We note that there is a one-to-one relationship between 
position and I in our structures, and so an enhancement 
factor of the form Fx{l) might not be as universal as 
Fig. ^ suggests. However, the energetic significance of 
the Laplacian terms, and the strong similarity between 
the form of the Laplacian and the LDA errors in Cxc (see 
Fig. 1^) , give us confidence in the physical importance of 
the Laplacian in describing inhomogeneity corrections to 
the LDA. Very recently, Becke Q and Perdew et al. ||] 
have constructed beyond-GGA functionals that contain 
the reduced Laplacian as a parameter. We consider this 
a step in the right direction and believe that the detailed 
local information on Uxc and Cxc in the strongly inhomo- 
geneous regime provided by our VMC calculations can 
be used to improve such functionals. 
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TABLE L Exchange-correlation energies (Hartrees per 
electron) and the LDA and GGA errors in this quantity for 
different values of the wave vector q. The statistical errors in 
E^J^'-^ are indicated. 



q/k% gJc^^^ AE^F^ AEg^^ 

1.11 -0.3289 ± 0.0001 -h0.0042 +0.0001 

1.55 -0.3127 ±0.0001 -0.0005 -0.0074 

2.17 -0.2882 i 0.0001 -O.OOfifi -0.0140 



FIG. 1. (color) The VMC and LDA n^c(r,r') for the 
q—1.55k% system plotted around an electron fixed at r (indi- 
cated by a white bullet) with r' ranging in a plane parallel 
to q: (top) electron on the slope; (bottom) electron at a den- 
sity minimum. The VMC hole is shown on the left and the 
LDA hole on the right. The charge density is also represented 
schematically. The color coding is only for rixc and varies be- 
tween -0.0500 (blue) and -0.0025 = (red). 
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FIG. 2. The upper graphs show e^f^ - e^^'^ along a di- 
rection parallel to q for two different strongly inhomogeneous 
systems. The lower graphs show the corresponding electron 
densities (light lines) and Laplacians (heavy lines). Distances 
are in units of the Fermi wavelength X% = 2n/k%. 




FIG. 3. The left panel plots the exact enhancement factor 

Fx against the reduced density gradient for three strongly in- 
homogeneous systems. The right panel shows the same quan- 
tity plotted against the reduced density Laplacian. 
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